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Abstract 

For the nonlinear matrix equations arising in the analysis of M/G/l-type and GI/M/l-type 
Markov chains, the minimal nonnegative solution G or R can be found by Newton-like methods. 
Recently a fast Newton’s iteration is proposed in : 14 ]. We apply the Newton-Shamanskii 
iteration to the equations. Starting with zero initial guess or some other suitable initial guess, 
the Newton-Shamanskii iteration provides a monotonically increasing sequence of nonnegative 
matrices converging to the minimal nonnegative solution. We use the technique in [?] to 
accelerate the Newton-Shamanskii iteration. Numerical examples illustrate the effectiveness 
of the Newton-Shamanskii iteration. 
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1 Introduction 

Some necessary notation for this article is as follows. For any matrix B = [b i;i ] € R nxra , B 
0 (B > 0 ) if bij > 0 (bij > 0 ) for all i,j; for any matrices A, B £ R nx ™, A > B (A > B) if a^ 
bij ( a,ij > bij ) for all i,j ; the vector with all entries one is denoted by e — i.e. e = (1,1, • • • , 1) T ; 
and the identity matrix is denoted by I. An M/G/l-type Markov Chain (MC) is defined by a 
transition probability matrix of the form 

Bo B\ B2 Bo 

C A\ A2 A3 

p — A 0 Ai A 2 

Ao Ai 

0 _ 

while the transition probability matrix of a GI/M/l-type MC is as follows 

0 

Ao 

Ai A 0 

A2 Ai 

where Bo € K m ° xm o and Ai € R mxm , respectively. N is the smallest index i such that A,;, for 
i > N, is (numerically) zero. The steady state probability vector of an M/G/l-type MC, if it exists, 


P = 


Bo C 
B\ Ai 
B2 A2 

Bo A3 
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can be expressed in terms of a matrix G that is the element-wise minimal nonnegative solution to 
the nonlinear matrix equation [7] 

N 

G = Y J A l G\ (1.1) 

i =0 

Similarly, for the GI/M/l-type MC a matrix R is of practical interest, which is the element-wise 
minimal nonnegative solution to the nonlinear matrix equation [5] 

N 

= ( 1 . 2 ) 

i=0 

It’s known that any M/G/l-type MC can be transformed into a GI/M/1- type MC and vice versa 
through either the Ramaswami HD or Bright [12] dual, and the G(R) matrix can be obtained 
directly in terms of the RiG) matrix of the dual chain. The drift of the chain is defined by 

P = P T P, (1-3) 

where p is the stationary probability vector of the irreducible stochastic matrix A = Aj, 

f3 = iAie. The MC is positive recurrent if p < 1, null recurrent if p = 1 and transient if p > 1 

— and throughout this article it is assumed that p ^ 1. 

Available algorithms for finding the minimal nonnegative solution to Eq. (11.11) include func¬ 
tional iterations [7], pointwise cyclic reduction (CR) [3], the invariant subspace approach (IS) [2], 
the Ramaswami reduction (RR) [4], and the Newton iteration (NI) [15] [6^ Ml I14j . For the detailed 
comparison of these algorithms, we refer the readers to [U] and the references therein. Recently, 
a fast Newton’s iteration is proposed in m and results in substantial improvement on CPU time 
compared with its predecessors. From numerical experience, the fast Newton’s iteration is a very 
competitive algorithm. 

In this paper, the Newton-Shamanskii iteration is applied to the Eq. m- Starting with a 
suitable initial guess, the sequence generated by the Newton-Shamanskii iteration is monotonically 
increasing and converges to the minimal nonnegative solution of Eq. (HU)- Similar with Newton’s 
iteration, equation involved in the Newton-Shamanskii iteration step is also a linear equation of the 
form BjXCj = E, which can be solved by a Schur-decomposition method. The Newton- 

Shamanskii iteration differs from Newton’s iteration as the Frechet derivative is not updated at 
each iteration, therefore the special coefficient matrix structure form can be reused. 

The paper is organized as follows. The Newton-Shamanskii iteration and its accelerated iter¬ 
ative procedure using a Schur-decomposition method are given in Section 2. Then M/G/l-type 
MCs with low-rank downward transitions and low-rank local and upward transitions are consid¬ 
ered in Section 3 and Section 4, respectively. Numerical results in Section 5 show that the fast 
Newton-Shamanskii iteration can be more efficient than the fast Newton’s iteration proposed in 
m- Final conclusions are presented in Section 6. 

2 Newton-Shamanskii Iteration 

In this section we present the Newton-Shamanskii iteration for the Eq. ED- First we rewrite 
(11.11) as 

N 

g(X) = J 2 A - xv ~ x = ° ( 2 . 1 .) 

v=0 

The function Q is a mapping from R mxm into itself and the Frechet derivative of Q at X is a linear 
map G x : R mxm —>. R m x m given by 

N v—1 

A(z) = EE A v X j ZX v ~ 1 - j - Z. (2.2) 

v=l j =0 
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The second derivative at A', Q ' x : K mxm —>. R mxm , is given by 

TV v—1 j— 1 TV v—2 v—2—j 

g'x(z 1 ,z 2 ) = EE^E X i Z 2 X j - 1 - i )Z 1 X v - 1 - j +EE^( E x^x”-*-*-*). 

v=2j=0 i= 0 v=2j—0 i— 0 

(2.3) 

For a given initial guess Gq,o, the Newton-Shamanskii iteration for the solution of 0(x) = 0 is 
as follows: 

for k = 0,1, • • • 

QGk o A ^’ s —1 — Q {Gk,s— l), G k ,s — G k: s—1 ~b Xk,s—li $ 1) 2, • , Tl }cj (2.4) 

Gfc+i = Gfc+ijO = Gk,n k ■ (2-5) 

Xfc jS _i is the solution to 

N u-1 N 

Afc )S _i E E A^ix^cr 1 -' 1 = E - G fc , s _ 1; 

i>—1 j=0 v=0 

which, after rearranging the terms, can be rewritten as 

TV—1 TV TV 

x M _r - e E ^ge^m-igi = e A - G is~i - Gfe,s-i- ( 2 . 6 ) 

j= 0 v=j+l n—O 

Following the notation of m, we define S k ,i = EE A j^t *> then the above equation is 

iV-1 JV 

(*S'fc,i — I)Xk, s -i + El = Gk , s -1 — El A ^k,s—ii (2-7) 

j=l u=0 

which is a linear equation of the same form EyE BjXC ■? = E as the Newton’s iteration step. It 
can be solved fast by applying a Schur decomposition on the matrix C, which is the m x m matrix 
Gk here, and then solving m linear systems with m unknowns and equations. For the detailed 
description for solving EyE BjXC 1 = E, we refer the reader to 13 s H4] , We stress that for 
Newton-Shamanskii iteration, the coefficient matrices are updated once after every rik iteration 
steps and the special coefficient structure can be reused, so the cost per iteration step is reduced 
significantly. 


3 The Case of Low-Rank Downward Transitions 


When the matrix Ao is of rank r, meaning it can be decomposed as Aq = AoF with Ao £ R mxr and 
r £ R rxm , we refer to the MC as having low-rank downward transitions. If Newton-Shamanskii 
iteration is applied to this case, all the matrices Xfc iS _i can be written as Afc iS _iF. This can be 
shown by make induction on the index s. Xo,o can be written as A'o^T and we assume that it is 
true for all X;j_i for l = 0 ,..., k and j = 1,..., s — 1. Hence G k ,s- 1 can be written as Gk, s - lL 

since G fcjS _i = EE E"=i x iJ -1 
m can be rewritten as 


E S— 1 y _ /•sr^k — 1 v 

j=i — vZ^z=o 2^j=i 


EE x k,j- i) r - Then 


TV TV 

Xfc )S _i = HoF + El A 3 G k,s-i G k,s-\F — G fe , s _ir + El A vG v k 1 X kt s-i 
i=i «=i 

N-l N 

+ E E ^G2- 1 -^ fc ,.- 1 Gjt- 1 G fc r > 

j =1 «=j+l 
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N 

= (i-X^gt 1 )" 1 

V = 1 

N N—l N 

x(A 0 + J2 A J G Z-I^k, s -1 - G fcjS _! ^ ^2 A v G°- 1 - j X kta _ 1 Gi- 1 G k ) r, 

j=i 3=1 «=i+i 

therefore Xfc jS _i can be decomposed as the product of an m x r matrix and an r x to matrix 

T. The inverse on the right-hand-side exists since 0 < J2v=i A vG v k ^ 1 < ’Y^^A V G V ~ 1 and the 
spectral radius of A V G V ~ 1 is strictly than one [3]. Therefore we will concentrate on finding 
as the solution to 

N N 

X k , s -i = A 0 + {^2 4? G i,s-i I)G k , s -i + ^ 2 AyG v k 1 X kiS -1 

3 = 1 v=l 

N — l N 

+ J2 J2 A v G°- 1 - j x kiB - 1 rG i k - 1 d k 

3= 1 v=j+ 1 

AT N- 1 

= Ao + (5Z ^3 G m-i - Skj+iXk,,-^TG k y, 

3 = 1 3=0 

which can be rewritten as 

JV-l N 

(S k ,i — I)X ktS -1 + ^2 Sk,j+iX kyS -i(TG k y = (I — ^ 2 AjG'l s 1 _i)G k ,s-i — Aq. (3.1) 

3 = 1 3 = 1 

We can use the Schur decomposition method in mm to solve the above equation. Different from 
the Newton’s iteration in [T4], the special coefficient structure can be reused here, thus saving the 
overall computational cost. We will repot the numerical performance of the Newton-Sliamanskii 
iteration in Section ?. 


4 The Case of Low-Rank Local and Upward Transitions 

In this section, the case of low-rank local and upward transitions is considered, where the m x m 
matrices {Ai , 1 < i < TV} can be decomposed as Aj = TAi with T £ M mxr an( j g To 

exploit low-rank local and upward transitions, we introduce the matrix U, which is the generator 
of the censored Markov chain on level i, starting from level i, before the first transition on level 
i — 1. The following equality holds based on a level crossing argument: 

N N 

U = Y J A l G i ~ 1 =J2M(I-U)- 1 A 0 y- 1 . ( 4 . 1 ) 

2=1 2=1 

For the case of low-rank local and upward transitions, we can rewrite U as 
N N 

u = J2 A *(( J - uy'Aoy- 1 = r[£ Adi - u^Aoy- 1 ] = r u, 

2=1 2=1 

which means U is of rank r, while G = (/ — U)~ 1 Aq is generally of rank m. 

Therefore we find U as the solution to 

N 

T(X) =X-J2 Mil - Xy'Ao)*- 1 = 0 , ( 4 . 2 ) 

2=1 
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and get G from G = (/ — £/) _1 Ao [TBJ 1111 . The Newton -Shamanskii iteration step for Eq. (14.21) 
is as follows: 

for k = 0,1, • • • 

•Fu k i = T(U k , a -i), T Y ks — i) s = 1,2, ■ ■ * , n k , 

Uk+l = Uk+ 1,0 = Uk,n k - 
Yfc jS _i is the solution to 

N 2—1 

Yk, a -1 - J2 Ai - U k )- l A 0 y-\l - U k )- x Y kia ^{{I - Uk^AoY^ 

i =2 j =1 

N 

= YtMV-^s-^AoY- 1 -U k , a -x. (4.3) 

2=1 

If we define = X^j+i A((/ — Uk) -1 Aq ) 1-1- ^(/ — U k ) _1 and rearrange the terms, Eq. (14.31) 
can be rewritten as 

JV-l N 

I'M-! ^ - U^AoY ~ Om-i) -1 ^) 4-1 “ 

j=i *=i 

which is of the form BjXC^ = E. This iteration enables us to exploit low-rank local and 

upward transitions. The iterates U k ,a = U k , s - 1 + Yjt )S _i, where solves Eq. m , can be 

rewritten as f7fc iS = TC/fc jS . This can be shown by make induction on the index s. It obviously 
holds for Uofl- Assuming U k , s -i = Tt/fc jS _i, from Eq. m we get 

N 2-1 

Y k , a -i = T[^A, ^((/ - U k )~ 1 A 0 ) i - 1 (I ~ U k )~ 1 Yk, a -i((I - U k )~ 1 A Q ) i ~ j 

2=2 j = 1 

N 

+ Y,MV - U^s.y^Aoy- 1 - u kta _ r], 

2=1 


which tell us that Yfc jS _i can be decomposed as TYfc )S _i, and the same holds for U k , s = U k , s -i + 
Yfc, s _i. Therefore from Eq. EM we will focus on finding Yk :S -1 as the solution to 


AT 2-1 

n, s _r - 

2=2 j= 1 

N 

= Y^Mv-^-^Aoy- 1 -u k , s -!. 

2=1 

Defining Ai((I — Uk) -1 Ao) 1-1- ^ (I — C4) _1 r, we can rewrite the above equation 

as 


N—l N 

Y k , a -1 - £fc,AM-i(( J - U^AoY = 55 A((I - - U k , a -u (4.4) 

i=i 


2=1 


which is of the form ]Tk =0 BjXC 3 = E. 


5 Convergence Analysis 

There is monotone convergence when the Newton-Shamanskii method is applied to the Eq. (D). 
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5.1 Preliminary 

Let us first recall that a real square matrix A is a Z-matrix if all its off-diagonal elements are 
nonpositive, and can be written as si — B with B > 0. Moreover, a .Z’-matrix A is called an 
M-matrix if s > p(B), where /?(•) is the spectral radius; it is a singular ill-matrix if s = p(B ), and 
a nonsingular M-matrix if s > p{B). The following result from Ref. [T7] is to be exploited. 

Lemma 5.1. For a Z-matrix A, the following statements are equivalent: 

(a) A is a nonsingular M-matrix ; 

(b) A- 1 > 0 ; 

(c) Av > 0 for some vector v > 0 ; 

(d) All eigenvalues of A have positive real parts. 

The following result is also well known mi. 

Lemma 5.2. Let A be a nonsingular M-matrix. If B > A is a Z-matrix, then B is a nonsingular 
M-matrix. Moreover, B < A -1 . 

The minimal nonnegative solution S for the Eq. (11.11) may also be recalled — cf. Ref. [tS; for 
details. 

Theorem 5.1. If the rate p defined by Eq. (El satisfies p ^ 1, then the matrix 

N v —1 

I - Y ^° V ~ 1 ~ 3 ) T ® A v° 3 

v=l j —0 


is a nonsingular M-matrix. 

5.2 Monotone convergence 

The following lemma displays the monotone convergence properties of the Newton iteration for the 

Eq. (11.11) . 

Lemma 5.3. Consider a matrix X such that 

(i) g(x) > o, 

(a) o <x <g , 

(Hi) I — Y^j=o(-^ V ~ 1 ~ 3 ) T ® AyXJ is a nonsingular M-matrix . 

Then the matrix 

Y = X- (Q' x )- l Q{X) (5.1) 

is well defined, and 

(a) G(Y) > 0 , 

(b) 0 <X <Y <G , 

(c) I — Li J2jZo(Y v ~ 1 ~ :i ) T ® AyYi is a nonsingular M-matrix . 
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Proof. g x is invertible and the matrix Y is well defined, from (iii) and Lemma ED Since 

N v—1 

[I-J2 ^2(X v - 1 - J ) T <g> AyXJ]- 1 > 0 

v=l j— 0 

from (iii) and Lemma 15.11 and Q(X) > 0, we get that vec(Y) > vec(X) and thus Y > X. From 
Eq. (15.11) and the Taylor formula, there exists a number 6, 0 < 9\ < 1, such that 

g(Y) = g(x) + g' x (Y-x) + ±g x (o 1 (Y-x),e 1 (Y-x )) 

= \g'x{0i{Y - X)MY - X)) 

> 0 , 


so (a) is proven, (b) may be proven as follows. From 

o = g[G) = g(X) + g x (G -x) + \y x {e 2 (g - x), e 2 (c - x)), 

where 0 < 0 2 < 1, we have 

~Gx(G-Y) = g' x (Y-X)-g' x (G-X) 

= -g(x)-g' x (G-x) 

= \g'x{02{G - X),6 2 {G - X)) 

> 0 , 

where the last inequality is from G — X > 0 by (ii). It is notable that 


(5.2) 


N v—1 

I-J2 51 (^"^'f ® A V X J 

v—1 j =0 

is a nonsingular M-matrix, so vec(G — Y) > 0 from Lemma IQI i.e. G — Y > 0. Now Y > X, 
so (b) follows. Next we prove (c). Since 0 < Y < G, 


N v—1 N v—1 

i - y, 5Z( y " _1 “ J ) T ® A - YJ > 1 - Y Y( GV ~ 1 ~ j ) T ® A -° 3 . 

v—1 j —0 v=l j —0 


and I — Xq= o( GV 1 J ) T ® A V G J ' is a nonsingular M-matrix. Consequently from Lemma 

15.21 I — J2v =i Y=o( YV ~ 1 ~ j ) T ® A v Y 3 is a nonsingular M-matrix. □ 

A generalization of Lemma 15.31 provides the theoretical basis for the monotone convergence of 
the Newton-Shamanskii method for the Eq. m- 

Lemma 5.4. Consider a matrix X such that 


(i) g(x) > o, 

(ii) 0 < X < G , 

(iii) I — Y=o(X v ~ 1 ~ j ) T ® A V X 3 is a nonsingular M-matrix . 

Then for any matrix Z where 0 < Z < X, the matrix 

Y = x- (g'z^gQ 0 (5.3) 


exists such that 
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(a) Q(Y) > 0 , 

(b) 0 <X <Y <G , 

(c) I — Ylv=l Y, V J Zo(y v ~ 1 ~ j ) T ® A v Yi is a nonsingular M-matrix . 

Proof. Since 0 < Z < X, 

N v —1 N v —1 

1-^2 J2(Z v ~ 1 ~ i ) T ® A v Z j > I — ^2 '52(X v ~ 1 ~ j ) T ® A v x j . 

v=l j —0 v—1 j —0 

From (iii) and Lemma l5.21 Q z is invertible and the matrix Y is well defined such that 0 < X < Y. 
Let 

Y = X- (Gx^OiX) , 

such that Y >Y from Lemma T5 . 2 1 As also Y < G from Lemma T5.31 (b) follows. Now 

/ - 51 ^2(Y v - 1 ~j) T <g> A V Y* 

v—1 j—0 

is a nonsingular M -matrix from Lemma [5..'il and Y > Y , therefore L— E^Li EjE) {Y v ~ l ~^) T ®A V Y^ 
is a nonsingular M-matrix from Lemma 15.21 Next we show (a) is true. From the Taylor formula, 
there exists two numbers 0 3 and 64 , where 0 < $ 3 , 64 < 1, such that 

G{Y) = g(x) + g' x {Y-x) + ±g'x{9 3 (Y-x),o 3 (Y-x)) 

= Q(X) + g' z (Y -X) + (G x - g' z )(Y -x) + ^g x (o 3 (v - x),e 3 (y - x)) 

= g" z {{Y - X), 04(X - Z)) + \g x {9 3 {Y - X),9 3 (Y - X)) 

> 0 , 


where the lat inequality holds since X — Z > 0 and Y — X > 0. □ 

The monotone convergence result for the Newton-Shamanskii method applied to the Eq. EH 
follows. 

Theorem 5.2. Suppose that a matrix Gq is such that 

(i) g(G 0 ) > 0 , 

fit) 0 < Go < G , 

(iii) I ~ E?= 1 E® A v Gi ts a nonsingular M-matrix. 

Then the Newton-Shamanskii algorithm {EZP ESP generates a sequence {Gk} such that Gk < 
Gk+i < G for all k > 0 , and linifc^oo Gk = G. 

Proof. The proof is by mathematical induction. From Lemma 15.41 

Gq = G 0 ,o < • • • < Gq :TIo = Gi < G , 
g(Gi) > 0 , 


N v—1 

1 ® a vG { 

v=l j—0 
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and 


is a nonsingular M-matrix. Assuming 


G(Gi) > 0 , 

Go — G 0 ,o < • • • < Go,n 0 = Gi < • • • < Gi- i , ni _ 1 = Gi < G , 


and that I — i 1 J ) T <S> A„A^ is a nonsingular M-matrix, from Lemma [Mil 


~iv-l-j\T 


G{G%+ 1 ) > 0 , 


Gi — Gifi < ■■■ < Gi tHi — G i+ i < G , 

and I — Ef=i Sj=o(^r+i 1-,J ) T ® A v G{ +1 is a nonsingular M-matrix. By induction, the sequence 
{Gfc} is therefore monotonically increasing and bounded above by G, and so has a limit G* such 
that G* < G. Letting i -» oo in Gj + i > G^i = Gi — (<7g.) -1 £(Gj) > 0, it follows that f7(G*) = 0. 
Consequently, G* = G since G* < G and G is the minimal nonnegative solution of Eq. OH . □ 


6 Numerical Experiments 

So, while more iterations will be needed than for Newton’s method, the overall cost of the fast 
Newton-Shamanskii iteration will be much less. 


7 Conclusions 
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